Effect of realistic astrophysical inputs on the phase and shape of the WIMP annual 

modulation signal 
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The orbit of the Earth about the Sun produces an annual modulation in the WIMP direct detection 
rate. If the local WIMP velocity distribution is isotropic then the modulation is roughly sinusoidal 
with maximum in June, however if the velocity distribution is anisotropic the phase and shape 
of the signal can change. Motivated by conflicting claims about the effect of uncertainties in the 
local velocity distribution on the interpretation of the DAMA annual modulation signal (and the 
possibility that the form of the modulation could be used to probe the structure of the Milky Way 
halo), we study the dependence of the annual modulation on various astrophysical inputs. We first 
examine the approximations used for the Earth's motion about the Sun and the Sun's velocity with 
respect to the Galactic rest frame. We find that overly simplistic assumptions lead to errors of up 
to ten days in the phase and up to tens of per-cent in the shape of the signal, even if the velocity 
distribution is isotropic. Crucially, if the components of the Earth's velocity perpendicular to the 
motion of the Sun are neglected, then the change in the phase which occurs for anisotropic velocity 
distributions is missed. We then examine how the annual modulation signal varies for physically and 
observationally well-motivated velocity distributions. We find that the phase of the signal changes 
by up to 20 days and the mean value and amplitude change by up to tens of per-cent. 
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I. INTRODUCTION 

Arguably the best motivated non-baryonic dark matter 
candidate is the neutralino (the lightest supersymmetric 
particle), and direct detection experiments are just reach- 
ing the sensitivity required to probe the relevant region 
of parameter space [1]. Since the expected event rates 
are so small ( O(10~^ — 10) counts kg~^day~^ see e.g. 
Refs. [2,3]) distinguishing a putative Weakly Interacting 
Massive Particle (WIMP) signal from backgrounds due 
to, for instance, neutrons from cosmic-ray induced muons 
or natural radioactivity, is crucial. The Earth's motion 
about the Sun provides two potential WIMP smoking 
guns: i) an annual modulation [4] and ii) a strong direc- 
tion dependence [5] of the event rate. In principle the 
dependence of the differential event rate on the atomic 
mass of the detector (see e.g. Refs. [2,3]) is a third 
possibility, however this would require good control of 
systematics for detectors composed of different materi- 
als. While direction sensitive detectors probably offer the 
best long-term prospects for the unambiguous detection 
of WIMPs, and development of such a detector is under- 
way [6], annual-modulation searches are already feasible 
using large detector masses [7] . 

If the local WIMP velocity distribution is isotropic 
then the annual modulation is roughly sinusoidal with 
a maximum in early June (when the component of the 
Earth's velocity in the direction of the Sun's motion is 
largest) and amplitude of order a few per-cent. The 
DAMA collaboration, using a detector consisting of Nal 
crystal scintillators, have detected an annual modulation 
with roughly these properties, which they interpret as 



a WIMP signal [7]. Assuming a standard halo model 
with a Maxwellian velocity distribution and circular ve- 
locity Vc = 220kms~^, they find a best fit WIMP mass 
= 52 GeV and cross-section C,cr-p = 7.2 x lO^^nb * 
with the 3-cr allowed region encompassing masses and 
cross-sections in the range 30 GeV < < 100 GeV ^ 
and 10"%b < CcTp < IQ-^nb [7]. 

Taken at face value this allowed region is incompatible 
with the exclusion limits from the Cryogenic Dark Mat- 
ter Search (CDMS) [9], Edelweiss [10] and Zephn I [11] 
experiments. While these exclusion limits depend rela- 
tively weakly on the WIMP velocity distribution (vary- 
ing by of order tens of per-cent for fixed Vc [12,13]) 
the annual modulation signal depends sensitively on the 
halo model assumed [14-21] and the region of WIMP 
mass-cross-section parameter space corresponding to the 
DAMA signal may be significantly enlarged if non- 
standard halo models are considered [14,15,19]. In par- 
ticular Belli et. al. [19] carried out an extended analy- 
sis of the DAMA data, for a large range of halo models 
and parameters, and found a significant enlargement of 
the allowed region to encompasses masses in the range 
30 GeV < < 270 GeV and cross-sections in the range 
10^^*^ nb < (ap < 6 X 10~^nb. However, as recently 



*Here (,~ Pxl (0-3 GeV cm~^) parameterizes the uncertainty 
in the local WIMP density, p^. 

^The lower limit va^ > 30 GeV is imposed by hand and 
represents the, somewhat model dependent, limit on the neu- 
tralino mass from accelerator sparticle searches (e.g. [8]). 
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pointed out by Copi and Krauss [20] and Forncngo and 
Scopel [21], the phase of the annual modulation can be 
significantly different if the WIMP velocity distribution 
is anisotropic and the time variation may not be close 
to sinusoidal [18,21]. In fact, in apparent contradiction 
to the results of Belli et. al. [19], Copi and Krauss [20] 
found, considering a similar range of halo models, that 
the DAMA annual modulation signal could not be made 
compatible with the exclusion limits from the CDMS and 
Edelweiss experiments. 

Motivated by this apparent conflict, and also the pos- 
sibility that if a WIMP annual modulation signal were 
detected its phase and shape [18,20,21] could allow us to 
probe the structure of the Milky Way (MW) halo, we 
study in detail the astrophysical uncertainties in the cal- 
culation of the annual modulation signal. In Section II 
we concentrate on the motion of the Earth relative to the 
Galactic rest frame. In Section III we discuss the local 
WIMP velocity distribution, in particular the importance 
of using physically and observationally reasonable mod- 
els. In Section IV we examine the effects of the modeling 
of the Earth's motion and the local velocity distribution 
on the phase and amplitude of the signal and the extent 
to which it can be approximated by a sinusoidal varia- 
tion, before concluding with discussion in Section V. 



Galilean transformation v ^ v = v + Vo(i), where Vc(t) 
is the Earth's velocity relative to the Galactic rest frame, 
which has two components: the Earth's orbit about the 
Sun and the Sun's motion with respect to the Galactic 
rest frame. 



A. Orbit about the Sun 

The Earth moves in a close to circular orbit, inclined 
at an angle of roughly 60° to the Galactic plane, with 
orbital speed = 29.79 km s^-'^. The annual modulation 
is mainly determined by the component of the Earth's 
motion in the direction of the Sun's orbit 



Vc, (t) = V0 + Vc sin 7 cos a{t) , 



(4) 



where a{t) = 2n{t - to)/T, T = 1 year, to ~ 153 days 
(June 2nd), 7 « 30° is the angle between the axis of 
the ecliptic and the galactic plane and t is measured in 
days. Since Vq ^ Vg, if the velocity distribution is close 
to isotropic, the differential event rate can be expanded 
in a Taylor series in cos a{t) so that, to first order, [4]: 



Ro{E) +RiiE) cos a{t) . 



(5) 



II. MOTION OF THE EARTH 

The differential WIMP elastic scattering rate due to 
scalar interactions is given by [2]: 



dR 
dE 



{E, t) = C^p 



P0.3 jmp + nij^ 



A^T{E, t)F^{q) 



(1) 



where po.3 = 0.3GeVcm~^ (and C is defined so that the 
local WIMP density is C, P0.3), is the WIMP scattering 
cross section on the proton, Vc is the local circular veloc- 
ity, A and F{q) are the mass number and form factor of 
the target nuclei respectively, E is the recoil energy of 
the detector nucleus, and T{E) is defined, so as to be 
dimensionless, as [2] 



T{E, t) = 



f 



(2) 



where /u(f) is the WIMP speed distribution in the rest 
frame of the detector, normalized to unity, and Umin is the 
minimum WIMP speed that can cause a recoil of energy 
E: 



\ 2m'^mA 



1/2 



(3) 



where ruA is the Atomic mass of the detector nuclei. 

The WIMP velocity distribution in the rest frame 
of the detector is found by making a, time dependent. 



with Ro{E) » Ri{E), i.e. the modulation is roughly 
sinusoidal. 

A commonly used [16,22,15,20] expression for the 
Earth's motion, in Galactic co-ordinates (X,Y,Z) where 
X is toward the Galactic center, Y in the direction of 
rotation and Z toward the north Galactic pole, is: 

Ve(t) = [— sina(t), cosa{t) sin 7, — cosa{t) cosj] . 

(6) 

This expression assumes that the Earth's orbit is circular 
and the axis of the ecliptic lies in the Y-Z plane [16]. 

A more accurate expression for the components of the 
Earth's velocity relative to the Sun can be found by us- 
ing the expression for the Sun's ecliptic longitude (from 
e.g. p77 of Ref. [23]), A = X 1°. 915 sing + 0.020° sin 2,g 
where L = 280°. 460 + 0°. 9856003 1 is the mean longi- 
tude of the sun, corrected for aberration, g = 357°. 528 -|- 
0°. 9856003 t is the mean anomaly (polar angle of orbit) 
and t is the time in days from 1200 Universal Time § on 
Dec 31st 1999, and transforming from ecliptic to Galactic 
co-ordinates (e.g. pl3 of Ref. [23]). Lewin and Smith [3] 
carried out his procedure and found: 



*For simplicity we ignore the Sun's motion with respect to 
the local standard of rest in this sub-section. 

^Universal Time is equivalent to Greenwich Mean Time to 
within an accuracy of seconds. 
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FIG. 1. The components of the orbital velocity of the Earth 
in Galactic X,Y,Z co-ordinates (ve,x top panel, Ve,Y middle 
panel and Ve, z bottom panel) found using the assumptions 
discussed in the text: i) assuming that the Earth's orbit is 
circular and the ecliptic lies in the X-Y plane (eq. (6), short 
dashed lines) ii) ignoring the ellipticity of the Earth's orbit 
and the non-uniform motion of the Sun in right ascension 
(eq. (8) and Ref. [18], long dashed) iii) including the ellipticity 
of the Earth's orbit but not the non-uniform motion of the 
Sun (Ref. [21], dot-dashed) and iv) including the ellipticity 
of the Earth's orbit and the non-uniform motion of the Sun 
(eq. (7) and Ref. [3], solid). Time is measured in days from 
noon on Dec 31st 2002. 

Ve(t) = We (A) [cos/3x sin (A - Ax), 

cos/3y sin(A — Ay), sin (A — A^)] , (7) 

where Wo(A) = [1 — esin (A — Ag)], e = 0.016722 and 
Aq = 13° ±1° are the ellipticity of the Earth's orbit 
and the ecliptic longitude of the orbit's minor axis re- 
spectively, and A = (-5°.5303,59°.575,29°.812) and 
A, = (266°.141,-13°.3485,179°.3212) are the ecliptic 
latitudes and longitudes of the of the (X,Y,Z) axes re- 
spectively. 

Fornengo and Scopcl [21] used a similar expression 
where the the non-uniform motion of the Sun in right 
ascension is neglected in the co-ordinate transformation 
so that the sin (A — A,) terms in eq. (7) above are re- 
placed by cos[uj{t — ti)] whore ti = (76.1,156.3,352.4) 
days and iv = 27r/(lyear). Gelmini and Gondolo [18] 
used a slightly simpler expression, which neglects the el- 
lipticity of the Earth's orbit and the non-uniform motion 
of the Sun in right ascension. In terms of the Sun's eclip- 
tic longitude, which under these approximations is given 
by X{t) = uj{t—ti) where ti is the spring (vernal) equinox, 
ti = 79.55 days **,: 



**This figure is for 2003 and is accurate to within about 0.04 
days. The time at which the Spring equinox occurs increases 
by 0.24 days a year (see e.g. Ref. [24]). 
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FIG. 2. The differences between the expressions for 
the velocity of the Earth in Galactic co-ordinates 
(A(ue, x), A(t;e,Y), A{ve,z), top, middle and bottom panel 
respectively) relative to the full expression including the el- 
lipticity of the Earth's orbit and the non-uniform motion of 
the Sun (eq. (7) and Ref. [21]). Line types as in Fig. 1. Note 
the larger scale used for A(t;e, y). 



ei sin X{t) — 62 cos X{t) 



(8) 



where ei and ei are in the Sun's direction at the Spring 
equinox and Summer solstice respectively: 



ei = (-0.0670, 0.4927, -0.8676) , 
62 = (-0.9931, -0.1170, 0.01032) . 



(9) 



These expressions for the components of the Earth's 
orbital velocity are displayed in Fig. 1 , and the deviations 
between the expressions in Fig. 2. In the plots time is 
measured in days from noon on Dec 31st 2002, but the 
correct time (measurcid in days from noon on Dec 31st 
1999) is used in the calculations. If time is erroneously 
measured from the beginning of the year rather than the 
beginning of 2000, then this leads to a horizontal shift 
in the curves of 0.27 days for every year elapsed since 
2000. The deviations in the Y component are reassuring 
small (of order a few per-cent), however assuming that 
the Earth's orbit is circular and the axis of the ecliptic 
lies in the Y-Z plane (eq. (6)) leads to large errors in the 
X and Z components. The other three expressions are in 
relatively good agreement for all three components. 



B. Motion of the Sun 

We now consider the motion of the Sun, which can 
be divided into two components: the motion of the lo- 
cal standard of rest (LSR) and the Sun's peculiar motion 
with respect to the LSR, VQ^pec- The "standard values" 
often used in WIMP event rate calculations date back to 
Kerr and Lynden-BeU's 1986 "Review of Galactic con- 
stants" [25], however since then the Hipparcos satellite 
has provided accurate measurements of the motions of 
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large numbers of nearby stars, allowing more accurate 
determinations of the Galactic constants, so the "stan- 
dard values" should be updated. 

The standard value used for the Sun's peculiar 
motion, with respect to the LSR, is U0,pec = 
(10, 15, 8) km s""'^ [25], while Lewin and Smith [3] use 
^^0,pec = (9, 12, 7) kms~^, with errors of order a few 
kms~^ in each component leading to an uncertainty of 
several days in the phase of the modulation [18]. The 
value determined more recently using data from Hip- 
parcos is W0,pec = (10.0, 5.2, 7.2) km s"^ [26,27], with 
the stated errors in each component being roughly 
0.5 kms~^. 

Assuming that the MW is axisymmctric, then the 
motion of the LSR is, by definition, {0,Vc{Ro),0) 
where Vc{Ro) is the circular velocity at the solar ra- 
dius. Kerr and Lynden-Bell combined a large number 
of independent determinations of the circular velocity 
and found Vc = 222.2 kms~^, with standard deviation 
20kms^"'^ [25]. The proper motions of Ccpheids mea- 
sured by Hipparcos allow accurate determinations of the 
Oort constants A and B [28], which lead to a value for 
the circular velocity, in terms of the Galactocentric dis- 
tance i?o: ^;c(-Ro) = (27.2 ± 0.9)(i?o/kpc)kms-i [27]. 
One possible caution is that this calculation assumes 
that stars move on circular orbits, which may lead to 
a systematic error as the MW is not exactly axisym- 
metric [27]. Accurate very long baseline interferome- 
try measurements of the proper motion of the SgrA* 
radio source at the Galactic center provides an al- 
ternative, marginally inconsistent, determination of Vc'. 
(30.1 ± 0.9)(i?o/kpc)kms-i [29,27]. In this case non- 
zero peculiar motion of SgrA* produces a systematic 
error. Using Rq = 8.0 ± 0.5 kpc from the most re- 
cent compilation of determinations of Rq [30], leads to 
Vc{Ro) = (218± 15) kms~^ from Cepheids and Vc{Ro) = 
(241 ± 17) kms-i from SgrA*. 

III. LOCAL WIMP VELOCITY DISTRIBUTION 

Data analyzes nearly always assume a standard halo 
model with an isotropic Maxwellian velocity distribution: 

/(v) = 3/2 3 exp (-vV<h) , (10) 

where Vc,h is the contribution of the halo to the lo- 
cal circular velocity, which corresponds to a spherically 
symmetric density distribution with p oc r"^ i.e. an 
isothermal sphere. Observations and numerical simula- 
tions indicate that galaxy halos are in fact triaxial and 
anisotropic however (see Ref. [13] for a review). As the 
size [14,16 21] and phase [20,21] of the annual modula- 
tion signal depend quite sensitively on the halo model 
assumed, the realistic modeling of the WIMP velocity 
distribution is crucial when extracting parameters or ex- 
clusion limits from data or comparing results from differ- 
ent experiments. 



The steady state phase-space distribution function of 
a collection of collisionless particles is given by the col- 
lisionless Boltzmann equation, and the velocity disper- 
sions, < V? >, of the system are calculated via the Jeans 
equations, which are found by taking moments of the col- 
lisionless Boltzmann equation (see e.g. [31]). Solutions 
to the Je an's equations ha ve the property that the ten- 
sor = (I'l — Vi){vj — Vj) is symmetric, so that at any 
point a set of orthogonal axes can by chosen such that c 
is diagonal (i.e. a?- = crf^Sij) [31]. However, as discussed 
in Refs. [31,32], there is no equation of state relating 
the components of the velocity dispersion to the den- 
sity, so solving the Jeans equations requires assumptions, 
which may or may not be physically reasonable, about 
the shape and/or orientation of cr. Evans et. al. [32] 
presented the logarithmic ellipsoidal model (which has 
potential (f>(x,y,z) = {v^/2)\n{x'^ +y'^p~'^ + z^q'"^)) 
which is the simplest triaxial, scale free generalization of 
the isothermal sphere. They argued that principle axes 
aligned with conical co-ordinates correspond to physical 
distribution functions and calculated the corresponding 
velocity dispersions. On the axes of the halo conical co- 
ordinates arc locally equivalent to cylindrical polar co- 
ordinates and the local velocity distribution can be ap- 
proximated by a multi-variate Gaussian with principle 
axes aligned with the (X,Y,Z) Galactic axes. This form 
has been widely used in calculations of the WIMP direc- 
tion detection rate [32,22,19,13,20,21], however a multi- 
variate Gaussian with arbitrary velocity dispersions need 
not correspond to a physically sensible halo model. For 
instance Evans ct. al. do not consider parameter values 
for which the ratio of any two of the velocity dispersions 
is more than 3:1, so as to avoid models which are afflicted 
by instabilities [32]. 

We will now briefly examine what sets of values of the 
velocity dispersions are reasonable for the logarithmic el- 
lipsoidal model, if the Sun is located on the intermediate 
axis of the halo (the deviations from an isotropic veloc- 
ity distribution are smaller on the major axis [32,13]). 
On the intermediate axis the velocity distribution can be 
written as 

1 / 1,2 vl vl \ 

^^^^ = (27r)3/2ara^az''''P " 2^ " 2^ J ' 

(11) 

with the velocity dispersions given by 

2 —4 
2 ^ <hP 

(2 + 7)(l-p-2 + (Z-2)' 
2 ^ <h(2g''-p'') 

2(1 - p-2 + ^-2) ' 

= 2(1-^,-2 +,-2) ' (12) 
where "p and q are constants which satisfy < < p < 1 
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and arc related to the axial ratios of the density distri- 
bution, /i^2, by 



2.5 



{p^q^ +p'^ - q^) 



q2 _|_p2 



■ p'^q^ 



-'2 — 



g (P g -p^ +q^ 

q,2 _|_ p2 _ p2q2 



(13) 



and 7 is a constant isotropy parameter, which in the 
spherical limit p = g = 1 is related to the anisotropy 
parameter 



/3 = 1 



<vi> 



2<v^> 



(14) 



by —7 — 213. Note that cr^ + cr^ = h ^'^'^ > a^. 

If we require Uj/3 < a\ < Scrj, so as to avoid models 
which might be affected by instabilities [32], then the 
velocity dispersions lie in the ranges 0.32 < a^^/Vc^h < 
0.71, 0.71 < cr0/uc,h < 0.95 and 0.24 < a,/vc,h < 2.'o. If 
we also require that < /? < 0.4, as found at the solar 
radius in simulated halos [33] , then Ur is further restricted 
to 0.71 < (Tr/wch < 0.92. In Fig. 3 we plot the values 
of CTr and 0-0 (cr^ ^ ^c.h ^ '^'^) which satisfy these two 
requirements. We still need to check that the axial ratios 
of the density distribution are in reasonable agreement 
with observed and simulated galaxy halos. There is no 
simple algebraic relation between the /j s and the ai s. 
If we require 0.6 < Ji < 1.0 and 0.3 < I2 < 1.0 (see e.g. 
Ref. [34] for a review of constraints on the shape of dark 
matter halos) then some of the previously allowed sets of 
velocity dispersions are now excluded. 

The logarithmic ellipsoidal model may not be a unique 
solution to the .leans equations (and there is no reason 
to expect the Sun to be located on one of the axes of 
the halo), but otherwise there is no reason for the princi- 
ple axes of the velocity distribution to correspond to the 
axes of the halo. In other words a multivariate Gaussian 
velocity distribution with axes corresponding to the axes 
of the halo and arbitrary velocity dispersions may not 
correspond to a physically reasonable halo model. While 
the constraints on the axial ratios and anisotropy of the 
halo we have imposed are not cast in stone, they illus- 
trate that only restricted sets of values of the velocity 
dispersions correspond to obscrvationally and physically 
reasonable halo models, and that the ratios of the veloc- 
ity dispersions can not be too large. Simulations provide 
further support for this argument. Helmi, White and 
Springel [35] examined the simulation particles within 
a 4 kpc box located 8 kpc from the center of a Milky 
Way like halo and found that, apart from a clump of 
fast moving particles from a late accreting subhalo, the 
velocity distribution was well approximated by a multi- 
variate Gaussian, with principal axis velocity dispersions 
in the ratio 1 : 1.08 : 1.27. If in fact the local dark mat- 
ter distribution is not smooth (see below for discussion 
of this possibility), then the local velocity distribution 
could not be approximated by a multivariate Gaussian 




FIG. 3. Physically and observationally reasonable values of 



<Tr and (ctz = {v, 



2 

c,h 



tIY^^). Inside the dotted lines the 



ratio of any two of the velocity dispersions is no greater than 
1:3 and inside the solid lines the anisotropy parameter /J is 
in the range < < 0.4. Restricting the ax;ial ratios of the 
density distribution would further rule out some sets of values 
inside the solid lines. 



with any set of velocity dispersions. We conclude that 
the large velocity dispersion ratios found by Forncngo 
and Scopel to produce extreme distortions of the annual 
modulation ((1:5: 4) and (10 : 1 : 3)) [21] are unUkely 
to correspond to a realistic halo model. 



It is possible that the local WIMP distribution may 
not be completely smooth, or in other words that the 
phase-space distribution has not reached a steady state. 
The density and velocity distributions of the "particles" 
in numerically simulated halos are relatively smooth at 
the solar radius [33,35], however WIMP direct detection 
experiments probe the dark matter distribution within 
a local sub-mpc region and even the highest resolution 
simulations have less than 100 "particles" per kpc^. Fur- 
thermore the first ncutralino clumps to form have mass 
more than ten orders of magnitude smaller than the 
smallest subhalos resolved in numerical simulations [36], 
and it is possible that the cores of some of these first 
small, high density, clumps could survive to the present 
day [36,13,37]. Stiff and Widrow [38] have recently de- 
veloped a technique which uses test particles to probe 
the velocity distribution at a single point within a nu- 
merical simulation. They find a local velocity distribu- 
tion consisting of a series of discrete peaks, and argue 
that the smooth background found previously in simula- 
tions may be, at least partly, a consequence of the use 
of a finite volume to measure the velocity distribution. 
Even if the local phase-space distribution bears no trace 
of the first WIMP clumps to form, more massive sub- 
halos which have been accreted onto the MW relatively 
recently, and have sufficiently eccentric orbits, could pro- 
duce coherent streams of high velocity WIMPs at the 
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solar radius [39,35]. 
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An additional complication is that the contribution of 
the visible components of the MW to the circular veloc- 
ity at the Solar radius is non- negligible, and may even 
dominate that of the halo (sec e.g. Rcfs. [40,33]). Self 
consistent halo models can in principle be built (taking 
into account various observations) using Eddington's for- 
mula [17,19], however building mass models of the Milky 
Way is a complex and ill-constrained process [40,27]; a 
small change in the observational constraints can lead 
to a large change in the properties of the dark halo. In 
Dchncn and Binncy's mass models, where the halo is de- 
scribed by a spheroidal density distribution, the contri- 
bution of the MW halo to the circular velocity at the solar 
radius lies in the range Vc,h{Ro) = 110 — 180kms~^ [40], 
while Moore et. al. find, for halos with a central density 
cusp p (X r~^-^, taking 190kms~^ < Vc{Ro) < 230kms~^ 
and using a smaller set of constraints, lOOkms"^ < 
^'c,h(^?o) < 120kms~^ [33]. We therefore conclude that 
a relatively conservative range of values is: 100 km s~^ < 
^^c, h(^o) < 200kms~^, and that the "standard" value, 
Vc,h{Ro) = 220kms~^, is probably too large. 

IV. RESULTS 

To compare theoretical predictions with experimental 
data requires the detector atomic mass, form factor and 
resolution, and the relationship between the nuclear re- 
coil energy and the energy detected, to be taken into 
account (see e.g. Ref. [3]). As the DAMA collaboration's 
raw data is not publicly available instead of working with 
the differential event rate we look at the detector inde- 
pendent quantity T{vniin, t), as defined in eq. (2). The 
relation between iimin and the detected energy is given, 
for various target nuclei and WIMP masses, in Fig. 4 of 
Ref. [21]. For reference, the DAMA experiment has a 
measured energy threshold of 2 keV which corresponds 
(for Iodine which dominates) to a recoil energy thresh- 
old of 22 keV, which for = 50 (200) GeV requires a 
minimum WIMP velocity Vmin = 309 (146) kms~^. 

A. Motion of the Earth 

We will start by examining the effects of the model- 
ing of the motion of the Earth on the phase and shape 
of the annual modulation signal. We use two fiducial 
halo models: the standard halo model (with an isotropic 
Maxwellian velocity distribution) and the logarithmic el- 
lipsoidal model [32] , with the Earth located on the inter- 
mediate axis and parameter values p = 0.9, q = 0.8 and 
7 = —0.62 corresponding to axial ratios 1 : 0.78 : 0.48 
and (3 = 0.4. The later model is intended as a specific, 
somewhat arbitrary but not unreasonable, example of an 
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FIG. 4. The dependence of r(vmin, t) on the approxima- 
tions used when calculating the Earth's orbit for, from top 
to bottom at t=0, Umin = 100,500,400,300 and 200 kms"^ 
with line typos as in Fig. 1, for the standard halo model (up- 
per panel) and for the fiducial triaxial model, see text for 
details, (bottom panel). Note the change in the phase of 
the annual modulation for Wmin = lOOkms"^ (see text and 
Ref. [41] for further discussion of the change in phase which 
occurs for small Vmin). We have fixed Vc,h = 150 kms"'^ and 
VQ,pec = (10.0,5.2, 7.2) km s"^ here. 

anisotropic halo model. In both cases we take Vc, h = 
150 km s^^. In Fig. 4 we plot T{vynin, t) produced by 
each of the expressions for the Earth's motion discussed 
in Sec. II, for five specific values of Wmin- As has long been 
known (see e.g. Ref. [41]), for small Wmin the phase of the 
modulation is reversed. Furthermore as fmin is increased 
the amplitude of the modulation reaches a maximum 
before decreasing again. For the standard halo model 
this maximum occurs at roughly 220 km s""'^, whereas for 
the fiducial triaxial model the maximum is at around 
250kms-i and r(200kms-\ t) ~ T(300kms-S t). In 



^^Our plots look superficially different to those of Forncngo 
and Scopel as they normalize to the event rate on January 
1st whereas we subtract the mean event rate. 
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FIG. 5. As Fig. 4 for the dependence of T(7;min, t) on the 
values used for the Sun's velocity with respect to the LSR 



in kms"^): (10.0, 5.2, 7.2) (Ref. [26,27], solid line), 



(0, 0, 0) (dotted line), (10, 15, 8) (Ref. [25], short dashed), 
(9, 12, 7) (Ref. [3], long dashed). The full expression for the 
Earth's orbit (Ref. [3] and eq. 7) is used here. 



Fig. 5 we plot T(wmin, t) produced by each of the sets 
of values for the Sun's motion with respect to the LSR 
discussed in Sec. II. 

To provide a quantitative comparison tables I and II 
contain the day at which T(?;iiiiii, t) (and hence the differ- 
ential event rate) is largest, tp, the maximum deviation 
from T(wmin, t) calculated using the full expression for 
the Earth's motion (cq. (7)), and the maximum devia- 
tion from a sinusoidally varying T{v^in, t) with the same 
mean, amplitude and tp, for the standard halo model and 
the fiducial triaxial model respectively. Tables III and 
IV contain the same data for r(wmin, t) found using each 
of the sets of values for the Sun's velocity with respect 
to the LSR. For reference the DAMA collaboration have 
carried out a fit to the modulation of their data using the 
function A cos [2-11 (t — tp)/T], and when T is fixed at one 
year they found tp = 144 ± 13 days [7]. 



Earth's orbit 


tpeak (days) 




*max 


zero-th order 


153 


-2.4% - -15.6% 


0.1% - -11.1% 


first order 


145 


2.6% - 19.0% 


0.3% - -6.6% 


GG 


148 - 149 


2.2% - 15.6% 


0.4% - -7.2% 


FS 


146 - 147 


2.2% - 15.8% 


0.5% - -6.7% 


LS 


147 - 148 


n/a 


0.7% - -5.1% 



TABLE I. The time, tp, at which T(t;inin, t) is largest, the 
largest percentage deviation from T(7;min! t) calculated using 
the full expression for the Earth's velocity, Aji^axi the largest 
percentage deviation from a sinusoidally varying T{vmin, t) 
with the same mean, amplitude and tp, A^^, for each of the 
expressions for the Earth's velocity discussed in Sec. II, for 
the standard isotropic halo model. Using only the compo- 
nent of the Earth's motion in the Galactic plane is denoted 
by "zero-th order", "first order" denotes assuming the axis 
of the ecliptic lies in the Y-Z plane (eq. 6), "GG" neglecting 
the ellipticity of the Earth's orbit and the non-uniform mo- 
tion of the Sun (Ref. [18] and eq. (8)), "FS" neglecting the 
non-uniform motion of the Sun (Ref. [21]) and "LS" including 
the ellipticity of the Earth's orbit and the non-uniform mo- 
tion of the Sun (Ref. [3] and eq. (7)). Where a range of values 
are given the quantity varies with Vmin and the left (right) 
hand value is for Vmin = 200 (500) krris^^. Here wo have fixed 
Vc,h = 150kms~^ and used VQ,pec = (10.0, 5.2, 7.2) kms~^. 



Earth's orbit 


tp (days) 


'-^max 


Asin 
^max 


zero-th order 


153 


-1.5% --11.6% 


0.1% - -4.7% 


first order 


147 - 142 


2.2% - 12.5% 


0.3% - -1.0% 


GG 


153 - 146 


2.0% - 9.6% 


0.4% - -2.3% 


FS 


151 - 144 


2.0% - 9.9% 


0.4% - -2.7% 


LS 


152 - 144 


n/a 


0.3% - -3.2% 



TABLE II. As Table I for the fiducial triaxial halo model. 



«0,pcc(kms ^) 


tp (days) 


A" 

'—'max 


^— *max 


(0, 0, 0) 


157 


3.1% - 19.3% 


0.7% - 5.6% 


(10, 15, 8) 


148 


4.0% - 30.0% 


0.6% - 4.9% 


(9, 12, 7) 


149 


2.7% - 20.0% 


0.6% - 4.9% 


(10.0, 5.2, 7.2) 


147 - 148 


n/a 


0.7% - 5.1% 



TABLE III. As Table I for the motion of the Sun relative 
to the LSR, with A^ax the maximum deviation from T(i;min) 
found using t;c3,poc = (10.0, 5.2, 7.2) km s^"*^, as found from 
Hipparcos data [26]. The full expression for the Earth's orbit 
(Ref. [3] and eq. 7) is used here. 
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W0,pcc(kms ^) 


tp (days) 


A hyp 
^— ^max 


A sinsoid 
^max 


(0, 0, 0) 


159 153 


2.4% - 13.8% 


0.1% - 2.6% 


(10, 15, 8) 


152 146 


2.8% - 20.1% 


0.2% - 2.5% 


(9, 12, 7) 


153 - 147 


1.9% - 13.6% 


0.2% - 2.8% 


(10.0, 5.2, 7.2) 


152 - 143 


n/a 


0.3 % - 3.2% 



TABLE IV. 
model. 



As Table III for the fiducial triaxial halo 



The time at which the event rate is maximum W varies 
by up to ten days depending on the expression used for 
the Earth's velocity, and depends on Wmin if the velocity 
distribution is anisotropic [20,21]. If only the component 
of the Earth's velocity in the Y direction is used, then the 
change in the phase which occurs for anisotropic halos is 
missed. The three sophisticated expressions produce re- 
sults which are in reasonably good agreement however, 
with maximum deviations of around a few per-cent for 
^'min ^ 400 km s^^. For larger Umin the exact form of the 
high energy tail of the speed distribution becomes impor- 
tant, so that using a different expression for the Earth's 
velocity produces a large fractional change in T(vmim t)- 
As the mean event rate is small for large Vmin, the abso- 
lute errors are small though and therefore not so impor- 
tant from an experimental point of view. 

The maximum deviations from a sinusoidal T{vynin, t) 
are no more than a few per-cent for Wmin ^ 400kms~^, 
as is expected since the error in neglecting the 2nd 
and higher order terms in the Taylor expansion is 
0{{vc.Y /vq)"^ ^ 0.01). The amplitude of the annual 
modulation is largest for large Vmm [4], as in the tail of 
the speed distribution the fraction of particles with speed 
greater than some fixed value changes significantly with 
Vf.(t), the Taylor expansion is then inappropriate and the 
deviation from sinusoidal variation becomes large (up to 
around 10%). For the fiducial anisotropic model we have 
chosen the deviations happen to be smaller than for the 
standard isotropic halo model, but this is not the case in 
general. 

The time at which the event rate is maximum is rel- 
atively weakly dependent on the Sun's velocity with re- 
spect to the LSR (provided that it isn't neglected en- 
tirely), however the maximum deviation from T{vrai-a-,t) 
found using the Hipparcos values for the Sun's velocity 
depends strongly on the component in the Y direction 
and can be large. Finally we note that combining a poor 
approximation for the Earth's orbit with an erroneous 
value for the Sun's velocity with respect to the LSR would 
produce even larger errors. 



I - 




**For simplicity we neglect the "flip" in the phase which 
occurs at small Wmin in this discussion. 



200 
t(days) 

FIG. 6. The dependence of T(umin, t) on the contribution 
of the halo to the local circulax velocity, Vc,h for the standard 
halo model for (from top to bottom) Vmin = 100, 200 and 
300 kms~^ and Vc,h = 100 (dotted line), 125 (short dashed), 
150 (solid), 175 (long dashed) and 200 kms"^ (dot-dashed). 

B. Local WIMP velocity distribution 

The region of WIMP mass-cross-section parameter 
space consistent with the DAMA annual modulation sig- 
nal depends strongly on the contribution of the halo to 
the circular velocity at the solar radius, i'c,h(-Ro) [14]. 
To explicitly illustrate the reason for this, in Fig. 6 we 
plot r(i;nnn, t) for i^niin = 100,200 and 300kms-i for 
values of Vc,-h in the range 100 — 200kms~^, for the stan- 
dard halo model. As Vc, h (and hence the typical WIMP 
velocity) is decreased the value of Wmin at which the mod- 
ulation flips phase and also that at which the amplitude 
of the modulation reaches a local maximum, are both 
smaller. For instance for v^.h = lOOkms"^ the local 
maximum occurs at roughly Wmin = 155 kms~^ compared 
with t;i„in = 220kms~^ for Uc, h = 150kms~^. 

As our aim is to examine how physical and observa- 
tional constraints restrict changes in the phase and shape 
of the annual modulation signal (rather than to carry out 
a detailed comparison with experimental data) we focus 
on the logarithmic ellipsoidal model, which reproduces 
some of the important features (triaxiality and veloc- 
ity anisotropy) of real galaxy halos. We examine the 
form of T{viam, t), with the Sun located on the interme- 
diate axis, for the sets of values of the velocity disper- 
sions that are unlikely to lead to instabilities and those 
which also correspond to halos with axis ratios and ve- 
locity anisotropy consistent with observed and simulated 
halos, as discussed in Sec. III. In Fig. 7 we plot the shift 
in the time at which T(wmin, is maximum, Atp, and 
the fractional change in the amplitude of the variation 
of T(vmin, t) relative to that for the standard isotropic 
halo model, for Wmin = 100,300 and 500 kms"-*^, fixing 
(y'ch = 150kms^^. In Fig. 8 we plot the fractional 
changes in the mean and amplitude of the variation of 
T{vmia, t), relative to the standard isotropic halo model. 
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FIG. 7. The shift in the phase, A(tp), and the fractional 
change in the amplitude of T(t;inin, t), A(amp), (relative to 
the standard isotropic halo model) for the logarithmic ellip- 
soidal model with the Sun located on the intermediate axis, 
for parameter values for which the ratio of any two of the ve- 
locity dispersions is no greater than 1:3 (dots) and for which 
also the anisotropy parameter is in the range < /3 < 0.4 and 
the axes ratios are 0.6 < /i < 1 and 0.3 < /2 < 1 (crosses). 
From top to bottom Wmin = 100, 300, 500 km s^^. The con- 
tribution of the halo to the local circular velocity is fixed at 
Vc,h = 150 kms"'^. 
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FIG. 8. As Fig. 7 for the fractional change in the mean, 
A(mean), and amplitude, A(amp), of T(t;inin, t) (relative to 
the standard isotropic halo model). 



If we consider all sets of velocity dispersions for which 

(Tj/3 < (Tj < 3(Tj, the shift in tp has magnitude of up 
to 40 days and is quite strongly dependent on Wmin- 
The change in the amplitude is of order 10s of per- 
cent for v^in = 100 and 300kms~^ and up to a fac- 
tor of 5 for Vaiin = 500 km s^^, while the change in the 
mean is even larger increasing from 10s of per-cent for 
Vmin = lOOkms"^ to morc than an order of magnitude 
at Umin = 500kms~^. The large change in the mean and 
amplitude for large Wmin occur because the logarithmic el- 
lipsoidal model has a more extended tail of high velocity 
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FIG. 9. As Fig. 7 for the logarithmic ellipsoidal model with 
t'c.h = 200 km s"'^, compared to the standard isotropic halo 
with Vc^h = 150kms~^. The diamonds denote the standard 
isotropic halo with Vc,h = 200kms~^. 



particles than the standard halo model [32,13] however, 
as noted in Sec. Ill, the mean event rate is very small for 
large Wmin- When only sets of velocity dispersions which 
correspond to halos with realistic axis ratios and veloc- 
ity anisotropy are considered the maximum shift in tp 
is roughly 50 % smaller and the the maximum change in 
the mean signal is roughly an order of magnitude smaller, 
which illustrates the importance of only considering ve- 
locity distributions that are physically and/or observa- 
tionally reasonable. Note that the standard isotropic halo 
model (which has /3 = and Ii = /2 = 1) lies at the edge 
of the observational constraints we impose. 

In Figs. 9 and 10 we plot the shifts in the mean, am- 
plitude and phase relative to the standard isotropic halo 
model with Vc^h = 150kms~^, for the same sets of ve- 
locity dispersion ratios but now with Vc,h = 200kms~^. 
The changes in the signal due to the change in v^. h arc of 
the same order of magnitude as the changes from varying 
the other parameters of the model. 

We do not find shifts in tp as large as those found by 
Copi and Krauss [20] and Fornengo and Scopel [21]. For 
many of the sets of parameter values we have consid- 
ered the shift in tp is large enough to be incompatible 
with the phase of the DAM A annual modulation signal, 
however. This suggests that if all three components of 
the Earth's velocity are included when calculating the 
event rate, then the consideration of anisotropic halo 
models may not lead to as large an increase in the re- 
gion of WIMP mass and cross-section parameter space 
corresponding to the DAMA annual modulation signal 
as found in Ref. [19]. As the DAMA data is not publicly 
available, however, it is not possible to confirm this. 

Finally we overview the effect of streams of particles 
on the annual modulation signal. A stream of particles 
with negligible velocity dispersion would produce a con- 
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FIG. 10. As Fig. 8 for the logarithmic ellipsoidal model 
with Vc,h = 200kms~^, compared to the standard isotropic 
halo with Vc^h = 150 km s"'^. The diamonds denote the stan- 
dard isotropic halo with Vc,h = 200 km s"'^. 



tribution to T(i>niin, t) which is a step the position 
and amplitude of which varies with time [39] : 



-^s (^min 5 t) 



2 CPO.3 Vs{t) 
0, for Vrnin > Vsit) 



, for Wmin < Vs{t) , 
t). 

(15) 



where ps is the density of the stream and Vs{t) is the 
speed of the stream in the rest frame of the detector, 
which is time dependent due to the time dependence of 
the Earth's velocity. The variation of the amplitude and 
position of the step depend on the relative orientation 
of the Earth's orbit and the orbit of the clump; for a 
clump of fixed velocity the closer the alignment of the 
clump's path with the Earth's orbit the larger the vari- 
ation in Vs{t) (in the, extremely improbably, case that 
the clump's path was perpendicular to the Earth's orbit 
there would be no variation). As outlined in Rcf. [39] 
if such a step were detected in the differential energy 
spectrum then the density, speed and direction of the 
stream responsible could, in principle (with a large num- 
ber of events and a detector with good energy resolution), 
be recovered. Multiple streams would produce multiple 
steps at different energies with different temporal vari- 
ations, which might be difficult to disentangle. Even if 
most of the dark matter is smoothly distributed a high 
velocity {vs ^ 500kms~^) stream of particles from a late 
accreting clump could produce a step in the differential 
event rate at large energy, the position of which would 
be modulated as discussed above [39]. However as the 
event rate decreases roughly exponentially with increas- 
ing energy this would only be detectable with very large 



''^The velocity dispersion would in fact be finite, but small, 
and would lead to a broadening in the step [39] . 



statistics (large target mass and exposure). A similar 
signal would be produced if there is a population of ex- 
tragalactic WIMPs [42]. 



V. DISCUSSION 

We have investigated the effect of uncertainties in as- 
trophysical inputs (the motion of the Earth with respect 
to the Galactic rest frame and the local WIMP veloc- 
ity distribution) on the calculation of the WIMP annual 
modulation signal. Accurate calculation of the shape and 
phase of the annual modulation signal requires all three 
components of the Earth's velocity with respect to the 
Sun to be taken into account properly. If the compo- 
nents perpendicular to the Sun's motion are neglected, 
the (energy dependent) shift in tp which occurs if the 
local velocity distribution is anisotropic is missed. Ne- 
glecting the motion of the Sun with respect to the Local 
Standard of rest, W0,poc, leads to an error of around 10 
days in tp, and uncertainties in VQ^pec lead to an error of 
order a few days in tp (as stated in Ref. [18]) and errors 
in the shape of the signal which grow from a few per- 
cent at low energies to more than ten per-cent at high 
energy. Approximating the modulation with a sinusoid 
with the same mean, amplitude and phase produces er- 
rors of up to ten per-cent. The errors in tp are crucial 
when comparing theoretical expectations with the mod- 
ulation observed by the DAMA collaboration [7] (which 
has phase tp = 144 ± 13 days). The errors in the shape of 
the signal induced by uncertainties in the Earth's motion 
are currently less important, but may become important 
for annual modulation searches with tonne size detectors 
(such as the planned GENIUS detector [43]), especially if 
we want to extract information about the local velocity 
distribution from an observed signal. 

If the local WIMP velocity is anisotropic then the 
phase [20,21], amplitude [16,15,18] and even shape [18,21] 
of the annual modulation signal can change. We have in- 
vestigated, focusing on the logarithmic ellipsoidal halo 
model, how the form of the annual modulation changes 
for parameter choices which correspond to physically and 
observationally reasonable halo models. We found that 
for reasonable sets of values for the velocity dispersions 
the shift in tp, which is energy dependent, can be up to 20 
days and the mean and amplitude of the signal change 
by tens of per-cent, at experimentally accessible ener- 
gies. It is possible that other halo models could produce 
larger changes in the annual modulation signal, however 
we have shown, in the context of this model, that re- 
stricting the choice of parameter values to those that 
are physically and observationally reasonable seriously 
restricts the changes in the signal. We have also argued 
that a multivariate Gaussian velocity distribution, with 
axes corresponding to the axes of the halo, and velocity 
dispersions with large ratios is unlikely to correspond to 
a physically reasonable halo model. 
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It is often assumed that the contribution of the lumi- 
nous components of the MW to the circular velocity at 
the solar radius is negligible so that Wc,h(^0) ^ Vc{Rq) ~ 
220 kms~^. This is not the case and in fact there is a large 
uncertainty in the value of Uc,h(^0) [40,33], which is im- 
portant since (as we saw in Sec. IIIB) the mean and am- 
plitude of the event rate (and hence the region of WIMP 
mass- cross-section parameters space corresponding to an 
observed modulation [14,19]) depend quite sensitively on 

Finally we discussed the possibility that the local 
WIMP distribution is not completely smooth [33,39,38]. 
If the local dark matter distribution consists of a num- 
ber of streams of particles with small velocity dispersion, 
then the event rate would contain a number of steps, 
the amplitude and position of which would vary with 
time [38]. Even if there is a smooth background WIMP 
distribution, high velocity streams from late accreting 
subhalos [39,35] may be detectable with good energy res- 
olution and a large number of events. 

In summary, it is important to properly take into ac- 
count astrophysical uncertainties (not just in the WIMP 
velocity distribution, but also in the motion of the detec- 
tor with respect to the Galactic rest frame) when calcu- 
lating the WIMP annual modulation signal. Analyzing 
data assuming a sinusoidal modulation with fixed phase 
could lead to erroneous constraints on, or best fit val- 
ues, for the WIMP mass or cross-section, even worse a 
WIMP signal could be overlooked. On the other hand 
using unrealistic halo models or parameter values could 
lead to overly restrictive exclusion limits or a mislead- 
ingly large range of allowed values of the WIMP mass 
and cross-section. 
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APPENDIX A: ERRATUM 

While the halo only contributes some fraction of the 
circular velocity at the solar radius, it is the total circular 
velocity Vc{Ro) = 220±20kms~^ (and not just the halo's 
contribution) which determines the velocity dispersions 
of the dark matter particles; the dark matter particles 
feel the gravity of all the components of the MW. Some 
of the values of vdRo) used for quantitative calculations 
in this paper are therefore too low. The main conclusions 
of this paper, regarding the effects of uncertainties in the 
modelling of the Earths's motion and the local WIMP 
velocity distribution on the phase and shape of the annual 
modulation signal, are unchanged however. 

I am grateful to Piero UUio and Joakim Edsjo for bring- 
ing this issue to my attention. 
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